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Abstract 

We apply path integration techniques to obtain option pricing with stochastic volatility 
using a generalized Black-Scholes equation known as the Merton and Garman equation. 
We numerically simulate the option prices using the technique of path integration. Using 
market data, we determine the parameters of the model. It is found that the market chooses 
a special class of models for which a more efficient algorithm, called the bisection method, is 
applicable. Using our simulated data, we generate some implied volatility curves. We also 
analyze and study in detail some of the characteristics of the volatility curves within the 
model. 
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1 Introduction 



Two popular instruments for option valuation in the market are the binomial option pricing 
and the Black-Scholes model The underlying stochastic nature of these models resembles 
quantum mechanical theories. In quantum mechanics, tools have been invented to solve and 
compute such stochastic quantities; one widely used instrument in this field is an intriguing 
integral called the path integration. The application of this method to option pricing has 
been mooted and elucidated in great detail in elsewhere Jl| . The underlying principle behind 
the method is essentially based on a generalized Black-Scholes model but unlike the Black- 
Scholes model, the new formalism can easily accommodate stochastic volatility and provide 
wider scope and greater flexibility for the investigation of market behavior. 

The most direct advantage of recasting the option pricing problem in terms of this Feyn- 
man path integral is that this formalism allows a new perspective for the trader and market 
analysts, and can lead to several new ways of computing exact, approximate and numerical 
solutions for the pricing of option. Stochastic volatility is naturally incorporated within the 
model. Stochastic volatility can introduce a high degree of nonlinearity within the option 
pricing problem. The path integral formalism can in principle handle such nonlinearity in 
an elegant manner. Moreover, there is a possibility that traders who are learnt this new tool 
can formulate new exotic options. 

There are many inherent similarities between the problem of pricing options in finance and 
solution of models in the physical sciences using stochastic approach. Indeed, this striking 
similarity have prompted Bouchaud and others |||,26] to apply successfully many mathe- 
matical tools previously used by physicists, like functional integration and scale invariance, 
to analyze problems in the financial markets. 

The celebrated Black-Scholes equation provides a simple and analytical formula for 
traders interested in plain vanilla European option. Equipped with the formula, analysis 
of the simplest option simply involves taking parameters estimated from historical data and 
working out the price of the option. Nowadays, even in the simplest scenario, market ana- 
lysts utilize option pricing using Black-Scholes formula with a twist. Instead of estimating 
the constant volatility required in Black-Scholes equation, which is often difficult and inac- 
curate, they usually compute the implied volatility which is necessary for the traded option 
price for a specific strike price to be consistent with the Black-Scholes formula. Such an 
analysis gives rise to a graph of implied volatility against strike price which normally ap- 
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pears as a smile or a frown. The decision regarding various investment portfolios depends 
largely on the difference between this implied volatility and the historical volatility Large 
implied volatility is taken by traders to mean that the option is overpriced. 

Numerous extensions on the Black-Scholes model have been proposed, see for instance 
Wilmott's book on financial derivatives [?] and references therein. Moreover, to explain the 
deviation of the call option from the Black-Scholes model, many refinements to the original 
model have been studied. For instance, one possible refinement is to relax the fixed volatility 
rate and replace it with stochastic volatility. Another refinement is to consider the stochastic 



model for the pricing of the security as a jump- diffusion process [pT^, p[| . Recently, it has 



been shown by Das and Sundaram [|TTJ that certain securities which follow a jump-diffusion 
process can exhibit stochastic volatility. There are essential differences between the last two 
examples of possible refinements. For the jump-diffusion case, the effect is enhanced within 
short terms whereas for stochastic volatility, the effect is more pronounced over a longer 
maturity time. Finally, one can also combine the two effects and consider stochastic interest 
rates. 

Constructing models for option pricing with stochastic volatility is an important issue 
since there are strong indications of stochastic variations in the underlying asset pricing and 
their derivatives from numerous empirical data. Several papers have recently attempted to 
study this pricing bias caused by stochastic volatility [?,p|,|TH, pj|, P^|. These studies have 
generally focused on numerical solutions of partial differential equations, Fourier inversion 
methods and the power series approximation techniques. The method of solving partial 
differential equations to study the option pricing has been considered to be the one of the 
general approach. However, a major distinct disadvantage of this technique is that most of 
the work done using this method is heavily computer intensive. The path integral formula- 
tion, on the other hand, offers an intermediate alternatives in many instances. It can include 
many of these techniques and offer new insights into the pricing of options. Moreover, the 
path integrals can yield results in a global approach involving the properties of the model at 
all times. 

In this paper, we investigate option valuation using the path integral approach. We gen- 
erate the volatility curves and study their behaviors. As the application of path integration 
to option pricing differs from the current theoretical method used for evaluating the option 
prices in the market, we shall briefly sketch the theoretical basis of our formulation in the 
next section. A detailed review and description of the path integral formalism in option 
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pricing will not be given here. A more elaborated account of the path integral formalism to 
option pricing can be found in the following reference PL 

In section |2], we discuss at length the theoretical framework of our model. The general 
algorithm using Monte Carlo techniques is given in the section [|. We perform the numerical 
simulations for our model in section |] and calibrate market data to the model in the same 
section. We find that the market seems to choose a special set of parameters. For this set 
of parameters, there exist more efficient algorithm called the bisection method. In section ||, 
we describe this algorithms and provide the necessary pseudocodes for the bisection method. 
We also briefly discuss the convergence of the program. 

Finally, in section [| we briefly discuss and study some characteristics of the implied 
volatility curves generated. Our simulation is based on Merton and Garman equation and 
consequently, there are several parameters which we can vary within the model. Besides the 
correlation coefficient between the security prices and the volatility, our model involves two 
other arbitrary constants; one related to the variance and the other related to the mean of 
the stochastic volatility. The initial stock price and the initial volatility can also be varied 
and studied. Using the path integral formalism, we generate computer data on option pricing 
by varying the various parameters in the model. Based on the option prices, we compute 
the implied volatilities and analyze the graphs of implied volatility against strike price. 



2 Theoretical Formulation 



Black and Scholes |4| laid the foundation for a quantitative analysis of European options. 
Since then, several extensions have been done and some of the original assumptions have 



been dropped. In an important paper, Merton |23[ ingeniously removed the assumption of 
constant interest rates and showed that an option can be priced in terms of a bond price. 
In the same paper, Merton also showed how the Black-Scholes formula can be extended to 
cover situations in which the volatility is a deterministic function of time. Indeed, based on 
a critical analysis of market data, Rubinstein has shown that the assumption of constant 
volatility is generally incorrect P5| . Research has also been done assuming different processes 
for the evolution of stock prices by Merton J24|, Cox and Ross |§ and Jones Cox and 
Ross U and Rubinstein |28| have solved the problem for the case when the volatility is a 
function of the underlying security price. 

Empirical evidence investigating the distribution of stock returns has shown mixed re- 
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suits. Kon [?] finds that the observed distributions are consistent with stochastic volatility 



while Scott BO] shows that the hypothesis that stock returns are distributed independently 



over time can be rejected. Bodurtha and Courtadon || and Hull and White |L5| also sup 



port the hypothesis of stochastic volatility. Considering these results, it seems reasonable to 
model volatility as another stochastic variable. 

2.1 The Stochastic Process With Volatility 

Several stochastic processes for the volatility have been considered by researchers. For ex- 
ample, Hull and White [i~4"|, Heston [[HJ and others have considered the process 



dV = (a + bV)dt + £V l/2 dz (1) 
where Q is white noise and V = a 2 . Baaquie M, Hull and White and others have considered 

dV = fiVdt + £Vdz (2) 

while Stein and Stein [?] consider 

da = -8{a - 6)dt + kdz (3) 

where 8 and 9 are constants representing the mean reversion strength and the mean value 
of the volatility respectively. We see that all the processes above except for (|3])0 follow the 
general form 

dV = (A + fiV)dt + iV a dz (4) 
The choice of A and \i is restricted by the condition that V > 0. 

2.2 The Merton-Garman Equation 

The process we are considering is 

dS = (pSdt + aSdzi (5) 
dV = (A + fiV)dt + iV a dz 2 (6) 



2 We can include this process if we add a term of the form jV 1 ^ 2 to the drift term 
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where 0, A, \x and £ are constants, V = a 2 and z\ and are Wiener processes with correlation 
— 1 < p < 1. Using Ito's lemma, we obtain the following expression for the process followed 
by a derivative fa dependent on the underlying security and the volatility of that security 

df * = (-m +<l)S ds + {x + flV) dv + — dsi- + pv u dsdv + —dv 2 ) dt 

We write it in this form to separate the stochastic and non-stochastic terms. 

We now consider two different options, f\ and / 2 on the same underlying security with 
strike prices and maturities given by K x , K 2 , Ti and T 2 respectively. We form a portfolio 

n = A + r 1 / 2 + r 2l s (8) 

so that 

du = (0i + r\e 2 + r 2 0,s)rft + (s : + r\s 2 + rvs)^i + (^i + v l m 2 )dz 2 (9) 

We have to get rid of the stochastic terms to ensure perfect hedging. Hence, we set 

Si + r\S 2 + T 2 aS = (10) 

^1+1^2=0 (11) 

to obtain 

r _ *i _ dh/dv 

1 v 2 df 2 /dV 1 1 
_y 1 df 2 dh _ dh/dv df 2 dfi 

2 ^ 2 dS OS df 2 /dV OS OS 1 ' 

Since the portfolio is now risk-less, it must increase at the risk-free interest rate by the 
principle of no arbitrage. In other words, we must have 

dli = rUdt (14) 
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Expanding II and simplifying, we obtain, after a separation of variables 
1 fdh q dfi VS 2 d 2 f\ 

dKJdv {-W + {x + » V) W + rS as + ~ as 2 " 
^ ? asdy 2 ov 2 J1 ) 
= WdV {-m +{x + » v) dv + r5 as + — as 2 " 

It is important to note that f3 is not a function of K±, K2, T\ or T 2 . This follows from the 
fact that the first expression is dependent only on K 1 and Ti while the second is dependent 
only on K 2 and T 2 . Hence, it is independent of all four variables. The term (3 is referred to 
as the market price of volatility risk. This is because the higher the value of /3, the more 
averse the investors are to take on the volatility risk. The reason this parameter is needed to 
price options with stochastic volatility and not for Black-Scholes pricing is that volatility is 
not traded in the market. Hence, it is not possible to perfectly hedge against the volatility 
even though it is possible to perfectly hedge against the underlying security price. Hence, 
investor risk preferences have to be considered when considering stochastic volatility or, in 
other words, risk-neutral valuation cannot be applied directly to volatility since volatility is 
not directly traded in the market. 

The parameter, f3, is difficult to estimate empirically and there is some evidence that it 
is non-zero [?]. To estimate this quantity, we consider the Cox, Ingersoll and Ross model 
where the consumption growth has constant correlation with the spot-asset return. This 
gives rise to a risk premium which is proportional to the volatility. We assume this model 
for simplicity as it has only the effect of redefining \i in the above equation. Henceforth, we 
shall assume that the market price of risk has been included in the Merton-Garman equation 
by redefining fi. Therefore, the Merton-Garman equation for the process we are considering 
is 

f + rS% + (A + ,V)% + \VS^ 2 + ff^S^L + ev^JL = rf (16) 
We introduce the variables S = e x and V = e y to simplify the calculations. In terms of these 
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variables, the Merton-Garman equation is 



,07 

dy 2 



( 17 ) 



For an option, we have f(T) = max(e x — K, 0), T being the maturity time. Hence, this is a 
final value problem. 

2.3 The "Straightforward" Solution when p = 

When p = 0, the solution for any volatility process, stochastic or non-stochastic is straight- 
forward. We make use of the theorem of Merton that the solution for a deterministic volatil- 
ity process is the Black-Scholes price with the volatility variable replaced by the average 
volatility. We can consider the stochastic volatility collection of a large number 

of deterministic volatility processes and the option price is then the average of the prices 
produced by each of the processes. In other words, if the volatility follows the generic process 
V(t) (where V may be stochastic), the option price will be given by 

poo 

C= / iSNid^V)) - Ke- rr N(d 2 (V))}V m (V)dV (18) 



where V m is the probability distribution function for the mean of the volatility ^ J^V(t)dt 
(which is a delta function for a deterministic process) and di(V) and ^(V") are the usual 

variate in Black-Scholes equation defined by dj = — — — - - — = - — -, j = 1,2, 

ay/T — t 

K is the strike price. Here, N(x) is the cumulative probability distribution function for a 
standardized normal variate. 

This intuitive result is derived in Scott . 



We will give two simple examples to illustrate this. First, let us consider a deterministic 
process. We will choose the process 

V = Voe*, < t < T (19) 



In this case, the probability distribution function of the mean of the volatility is given by 



e^ T - 1 ' 

V m = 6[V-V =-) (20) 
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giving us the Black-Scholes result with a replaced by yVo e ^ T 1 • 

For a stochastic volatility process, we choose^A = p = a = in eq(||) to obtain 

dV = £dz, 1/(0) = V , < t < T (21) 

where Q represents white noise. The distribution of the mean of V during the time interval 
(0, T) is given by 

V m ~N (Va, Sp\ (22) 
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Hence, the option price is given by 



C = ^2^r [SN ( dl(y)) ~ Ke ~ rTN ^ V ^ ex P (- 3( ^T° )2 ) dV - (23) 
2.4 An Extension to Merton's Theorem 

The case p = corresponds to Merton's Theorem. We now extend Merton's theorem to the 
case of non-zero correlation for the stochastic process of the volatility that we are investi- 
gating. The present value of the option is given by 

/"OO 

f(x, y,T) — / (x, y | e~ kr \ x')f(x', T)dx' (24) 



and (x,y \ e Ht \ x') is given by equation (|90|) with 5o and given by ( p3|) and ([54]) respec- 
tively. Now, since J DYe So describes the probability of a specific path for y (we show this 



in detail in chapter 6), we see that the propagator can now be written as in equation ( [59]) 
with uj, 9, 7] and £ being functionals of this path and v being the final value of V 3 ^ 2 ~ a for 
the path. Hence, if we generate paths for y according to @, the option price is given by 

/ poo Si(x,x',U!,9,ij, C, v) \ 

M= \L" im^ ( "'~ K)+ ) (25) 

(since the payoff of the option is given by (e x> — K) + ) with the average taken over the paths 
for V. Since the propagator is in the form of a Gaussian, we can perform the integration 
over x 1 to obtain 

/(*, y,t)=( SN( Sl ) exp {- P \ - (V^ ~ Vf-) -(fe-f V 



+ PlA-Ke-^N(s 2 ) 



3 This is not a realistic process as P(V < 0) > while V is obviously non-negative. However, it might be 
a reasonable approximation for relatively short times for which P(V < 0) is negligible. 
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where Vi and Vf are the initial and final volatilities of the path respectively and si and s 2 
are given by 

In (|) + rr + |(1 - 2p> + ^ (tf »"« - ^) + *0 + f rj - f ( ^ 
s i — / (27) 

S2 =Sl _ ^(1-^2)^ (28) 

It is easy to verify that equation (^(j) is the same as the Black-Scholes equation for any single 
volatility path when p = 0. 

When a = |, several simplifications occur. We see that 9 = ( = r are known and r] — u. 
In that case, ( p6| ) reduces to 

5iV( Sl ) exp (- f £ + Ejfju - | (V, - V» - p (^^-) r) - Ke^N(a 2 )^ (29) 
where si and S2 now have the relatively simple forms 

■n(#) + (r + »( 3 ^)) T -(^-f-Q-'+f( y '- y /) 

S2 =Sl _ y/^-ffuj (31) 

Hence, we see that we have a straightforward solution for a — | even when the correlation 
is not zero. 

When « = | and A = 0, we obtain a similar simplification since ( = to and rj = r. In this 
case, we obtain the following expression for the option price 

f(x, y, t) = (SN{ 81 ) exp ({p - flu + 2 In Q|) + 2 y ) ) ~ ^ e ~ rrAr (^)) (32) 
and si and s 2 are now given by 

+ (r + f)r+(-p 2 -f + §)u, + flnfe) 
si = . — - (33) 

s 2 = Sl - ^(I-P 2 )^ (34) 

For the case considered in Baaquie ffl, we have A = and a = 1. In this case, we have 
rj = C — Jq e y / 2 dt which gives us 

/(*, V, t) = (SNM exp (y^ - V// 2 ) -f(,-f),)- ^e-iV( S2 ; 

(35) 
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where Si and S2 are now given by 



7 \ /' - T J '/ 



In (#) + rr + ( V + |) a, + f (tf" - + £ (, 
si = . — (36) 

S2 = 8l - y/(l-p 2 )LU (37) 

which is somewhat more complicated since two functionals, uj and 77 of the volatility path are 
involved. In this case, however, a perturbation analysis can be used to derive an approximate 
form for the probability distribution functions of the functionals. Due to this fortunate 
occurrence, a series solution to this problem can be obtained. 

The probability density function (pdf) for the functionals is a very difficult quantity to 
obtain. The probability density function for JJ" Vdt was obtained for the special case a = 1/2 
in Stein and Stein [?][]. Stein and Stein [?] have used this probability distribution function 
and the "straightforward" solution for p = to get an analytic form of the solution for this 
case. We now see that the result can be extended to non-zero p if we can find the joint 
probability density function of this functional and Vi — Vf. While the individual probability 
distribution functions can be obtained (the pdf for uj is obtained in Stein [?] and the pdf for 
Vi — Vf is trivial), they are not independent. 



2.5 Risk-Neutrality 

We show that the expected value of the underlying security S whose initial value is So is 
given by Soe rt after time t has elapsed for a large class of stochastic processes including the 
one we are considering in this thesis. In other words, we show that A = e~~ rt S is a martingale. 

We first change variables from S to A in (|5|) (changing to r in accordance with risk- 
neutral valuation) to obtain 

dA = Ae y/2 d Zl (38) 

(where z\ is the time integral of W and hence a Wiener process) where y may depend on A 
(y can be stochastic). We now consider the more general process 

dA = f{A,y)dzi (39) 

The original solution for simple Brownian motion was obtained way back in 1944 by Cameron and 
Martin @ 
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We note that E[dA] = so that E[A(t + dt)\A(t) = A ] = A . In other words, we see that A 
is a martingale. Hence, we have shown the result. In general, a martingale process cannot 
have a drift term. 

While the result is simple, it has important consequences. We note that risk-neutrality 
alone cannot determine any constraints for the volatility process. Any volatility process 
whatsoever satisfies risk-neutrality. 



3 Numerical Algorithm 

3.1 A Short, Quick Reminder of the Monte Carlo Method 

The Monte Carlo algorithm to integrate 

/ f(X)dX (40) 

J A 

where X can be (and, in fact, usually is) a multi-dimensional variable and A is a subset of 
the domain of X requires us to split f(X) = g(X)p(X) so that J A p(X)dX = 1 (in other 
words, so that p(X) is a valid probability density function). The algorithm then states that 
an estimate for the integral is given by (g(Xi)) = 4 Yli=i 9{Xi) where the configurations Xi 
are generated randomly according to the probability density function p(X). 

The error of the Monte Carlo method goes as ^= as a consequence of the central limit 
theorem as long as there is no correlation between the configurations produced. (Though 
in general this condition is difficult to satisfy, we shall see later that we can easily satisfy it 
for this case). While this error may not look very impressive, it is often the best that can 
be managed for X which have a large number of dimensions. For the present problem, we 
have a very large number of dimensions (in fact the exact problem has an infinite number of 
dimensions) and the Monte-Carlo method is the most practical one available for it. 



3.2 A Monte Carlo Method for this Problem 



The Monte Carlo based numerical approaches of Hull and White [14 , 16 ] , Finucane [ 12 1 and 



n 



N ' 



Mills [25] are all generalized forms of the Binomial tree in which for discrete time t n 
the stock price x n and volatility y n are considered as random variables and both of which 
are updated. For the case of iV-steps, we need to update N 2 variables for obtaining a new 
configuration. In eq(|90"D, we have a drastic simplification since the path integral over the Xi 
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variables for the general case of p 7^ has been evaluated exactly by analytical means. Hence 
in basing numerical simulation on eq(pU|), one needs to generate new configurations for only 
the y n random variables, namely iV variables, reducing the computational time required by 
a factor of N. 

For this problem, we choose the following probability density function 

= n -7^ e ^ 



where Y is the set of variables yi (and is hence iV-dimensional) and So is given in (|93|). 
Hence, we see that 



e Si 



9(Y) = . — (42) 

/ 27re(l - p 2 ) Y:f =1 z yi 



where S\ is given in (]83Q since the integral we are performing is 



'N-l 



/ DY TT e^ 1 "^ (43) 



. i=i 
where 

JV-l 



7T6 



We now have to produce configurations Y with the probability distribution p(Y). While 
p(Y) looks rather complicated, it has a simple interpretation. It is the probability distribu- 
tion for a discretized random walk performed by y. To see this, let us first use Ito's lemma 
to find the process followed by y. We find, from eq(proc2) and eq(variable), 

/ t2 2y(a-l)\ 

dy = \ e - y + 2 dt + £e v{a - l) Qdt (44) 



2 

We can now discretize the process using Euler's method to obtain 

/ t2 2 Vl (a-l)\ 

8 Vi = f Ae~ w + /i - j e + fe^-^ZVe (45) 

where <5?/j = yi — yi-i, e is the time step and Z is a standard normal variable. Since we are 
using r = T — t as the time variable, the time step is actually — e. Hence, 5yi is a normal 
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random variable with mean (— Ae Vi — fi + £ 2 e 2yi ( a 1 - ) /2) e and variance £ 2 e 2yi( - a x >e and the 
probability density function is given by 

/i=?r^=exp [-^- d -^+Xe-^ + n-^- I (46) 



Hence, the joint probability density function for the discretized process is given by 

f=UA= nwk° (47) 



4 = 1 \ 1=1 v ^ 

which is the same as ((4~l|). 

In this simulation, we will use Euler's method to find the volatility paths since these 
paths are generated with the requisite probability distribution. 

The algorithm to find a Monte Carlo estimate of the propagator p = (x,y \ e~ Hr \ x') is 
as follows 

1. p := (Initialization) 

2. For i := 1 to N 

3. Generate a path Y for y using (fHj) 

4. p ■= p + g(Y)/N (where g(Y) is defined in (PD) 

5. End For 

The paths must be generated backwards starting from y^ which is the initial value of In V 
to obtain all the y^. Since the equations are time symmetric (after, of course, reversing the 
drift terms), this presents no problem. This will have to be repeated for all the x' that we 
wish to integrate over. However, during implementation it is found to be more advantageous 
to generate the paths only once, storing the important terms t = Ne, t\ = J2iLi eVl anc ^ 

f2 ^f e ^,(^ +Ae -», + ,_^-n 

i=l ^ ' 

which are sufficient to determine S\ once x' is given (Si = — 2 e(i- P 2 )t 1 { x — x' + (r — q)t — e (y + t 2 ))' 
where q is the annualized dividend. That Si can be computed using this limited information 
is fortunate as storing all the paths explicitly would require a very large memory (10MB for 
10,000 configurations as compared to 160kB when storing only the essential combinations of 
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terms). The alternative of generating paths for each value of x' (as in the naive algorithm 
above) is unacceptable due to the very large run time required using this approach. 

The propagator must finally be integrated over the final wave function to obtain the 
option price. The accuracy of this numerical quadrature depends on the spacing h between 
successive values of x' . This means that we have to find the propagator for several values 
of x' to obtain reasonable accuracy which is computationally very expensive. We found it 
better to find the propagator using the above Monte Carlo method for only about 100 equally 
spaced values of x' over the range of the quadrature and using cubic splines to interpolate 
it at the other quadrature points. This produces excellent results (as we shall see later) as 
the propagator seems to be an extremely smooth function of x' . 

Hence, the revised algorithm is of the form 

1. For i := 1 to N 

2. Generate a path Y for y using 

3. Store t\ and £ 2 for the path 

4. End For 

5. For x' := beginning of range to end of range 

6. Find the Monte Carlo estimate for the propagator at large intervals of x' using t\ and 
i 2 from the paths. 

7. End For 

8. For x' := beginning of range to end of range 

9. Find the propagator at small intervals of x' using cubic spline interpolation over the 
values of the propagator found previously and integrate over the final wave function 
(payoff) 

10. End For 

11. Return option price 

The numerical quadrature was performed using Simpson's rule. While Simpson's rule is a 
relatively low order method, it was deemed sufficient as the error in quadrature is negligible 



43) 
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compared to the Monte Carlo error in the propagator. For the special case of the European 
call option, it is possible to integrate analytically over the splines as the payoff is piecewise 
linear. However, this was not done as firstly, the error in the quadrature was negligible 
compared to the Monte Carlo error and secondly (and more importantly) we wanted to 
make the program general so that it is able to price any derivative based on the underlying 
security which has only one payoff date. 

For the special case of functions which can be expressed as piecewise polynomial functions, 
the integration can be carried out analytically and we can then just generate the paths and 
find the expectation over the generalized Black-Scholes prices. While this method has the 
advantages of being elegant, easy to program and completely eliminating the error due to the 
integration, it has the disadvantage that we have to carry out the integration first before we 
can write the program which will be specific to the given payoff. We have implemented this 
approach for call options but found no significant improvement in run time (it runs about 
1% faster). 

The above algorithm may appear unstable in the limit e — > since in this case — = 
O i^\J~^\- Fortunately, in the case of our simulations this is not a problem as our step size 
was always large enough to avoid this problem. If this is a problem in any simulation, this 
can be easily handled as the propagator can always be written as a functional of the paths 
as given in (|25|). We can then generate paths for V or y using (|6|), store the functionals 
(we will now need 5 terms t\ = uj (which is the same as above), t<2 = 8 = Yli=i e Vi<yl ^ 2 ~ a \ 
t 3 =r] = E 2 =ie yi(3/2_Q) , U = ( = J2?=i eVAa ~ 1/2) and *5 = e yo( - 3 / 2 - a >). We can then proceed 
to find the propagator as above (of course, now writing Si as a function of ij, i = 1, . . . , 5) 
and integrate over the final payoff. For the special case of the European call option, one 
can make this method slightly more efficient by averaging over the generalized Black-Scholes 
prices ( p6|) . We did not use this method since it takes longer to compute the propagator 
using five terms and because the memory requirements are higher. 

4 Numerical Results for the General Case 

We performed simulations on 90 day options setting an initial volatility of 25% per annum 
which is a figure comparable to the historical data. We performed simulations for a = 0, |, | 
and 1 with correlations ranging from -0.5 to 0.5 in steps of 0.1. We set A = fi = for most of 
the simulations. Some simulations with mean-reverting volatility processes were performed 
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for the purposes of investigating the effect of mean-reversion on option prices. 

The following parameters have the following values unless otherwise stated : Sq = 
100, r = 5%, q = 0%, t — 90 days, A = and p = 0. Most of the simulations have been per- 
formed using 128 time steps and 500,000 Monte Carlo configurations. The exceptions have 
been for a = | where we have used 512 time steps and the simulations we have performed 
for the other values of a to check that the number of time steps is sufficient. The error bars 
in all the graphs refer to the standard error obtained from the Monte Carlo simulation. 



4.1 The Effect of p on Option Prices 



Implied Volatility Curves Implied Volatility Curves 



for a = 0.0, \ = 0.02, o = 0.25 for a = 0.0, \ = 0.02, o„ = 0.25 



E 




$80 $90 $100 $110 $120 $130 "' $96 $98 $100 $102 $104 

Strike Price (Initial Security price = $1 00) Strike Price (Initial Security price = $1 00) 



Figure 1: Implied volatility curves showing the effect of p on option prices when a = 0. We 
can see that positive p leads to an increase in the option price when the strike price is high 
and a decrease when the strike price is low while negative p has the opposite effect. 

The correlation p has a very large impact on the implied volatility curve irrespective of 
the values of the other parameters. When p = 0, we see that the implied volatility is in the 
form of a smile with the minimum near the present value of the underlying security (Sq). As 
p increases, we find that the implied volatility increases for large strike prices and decreases 
for small strike prices. This can be easily seen in the graphs in figures [I], |2|, |3] and ^. We can 
also see that the deviation of the implied volatility from the initial volatility is much higher 
when the correlation is non-zero. These results are consistent with those reported in Hull 
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Implied Volatility Curves Implied Volatility Curves 



for a = 0.5, % = 0.1, o = 0.25 for a = 0.5, \ = 0.1, o = 0.25 




$80 $90 $100 $110 $120 $130 "' $96 $98 $100 $102 $104 

Strike Price (Initial Security price = $1 00) Strike Price (Initial Security price = $1 00) 

Figure 2: Implied volatility curves showing the effect of p on option prices when a — \. We 
can see that positive p leads to an increase in the option price when the strike price is high 
and a decrease when the strike price is low while negative p has the opposite effect. 



and White [15[], Heston ||13|| , Johnson and Shanno and Scott |30| . 

This can be explained in terms of the propagator obtained by the Monte Carlo sim- 
ulation. From figure [|, we see that the propagator for positive correlation is greater for 
very large x' — x (where x = In So and x' = In S) as compared to the propagator for zero 
correlation. This essentially means that the probability for the underlying security price 
reaching very large values is higher when the correlation is positive. Hence, the price of 
options with a very large strike price is larger when the correlation between the processes for 
the underlying security price and the volatility process is positive. For negative correlation, 
we see that the propagator for small, positive x' — x is greater than that for zero or positive 
correlation. Hence, for options whose strike price is smaller, there are two competing factors, 
the propagator for x' slightly larger than x and the propagator when x' » x. For relatively 
large strike prices, the latter is more important and the implied volatility is higher when the 
correlation is positive while for relatively small strike prices, the former is more important 
and the price when the correlation is positive is lesser than for zero or negative correlation. 
The same analysis can be performed for the case of negative correlation and is consistent 
with the simulated results. 
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We can also intuitively examine why the propagator has this form when the two processes 
are correlated. If the two processes are positively correlated, we have two possibilities. If the 
price initially increases, so will the volatility and hence, large price increases become more 
likely while small price increases become less likely. On the other hand, if the price initially 
decreases, so will the volatility and the price is more likely to remain at that value. Hence, 
for positive correlation, the propagator must be higher for prices slightly lower and much 
greater than the initial price and lower for prices slightly larger than the initial price. The 
reverse holds true for negative correlation. This naive reasoning is fully borne out in figure 



We were unable to simulate frowns or any sort of kinks in the implied volatility curve 
for any values of the parameters. This seems to suggest that stochastic volatility even with 
arbitrary correlation places some constraints on the shape of the implied volatility curve. We 
can use this to check whether the hypothesis that the volatility is stochastic is reasonable or 
otherwise. We note that the empirical implied volatility curve used for our calibration does 
satisfy this criterion. 

Implied Volatility Curves Implied Volatility Curves 



for a = 0.75, ^ = 0.2 p = -0.9 to 0.9 



0.255 



for a = 0.75, \ = 0.2 p = -0.9 to 0.9 




Figure 3: Implied volatility curves showing the effect of p on option prices for a — 1. We 
can see that positive p leads to an increase in the option price when the strike price is high 
and a decrease when the strike price is low while negative p has the opposite effect. 



19 



0.280 



0.260 



0.240 



0.220 



Implied Volatility Curves 

for a = 0.75, (, = 0.2 p = -0.9 to 0.9 
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Implied Volatility Curves 

for a = 0.75, \ = 0.2 p = -0.9 to 0.9 






98.0 
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Strike Price 



102.0 
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Figure 4: Implied volatility curves showing the effect of p on option prices for a = 1. The 
curves for the different values of p are in ascending order according to the slope (in other 
words, the slope increases monotonically with p). Hence, we can see that positive p leads to 
an increase in the option price when the strike price is high and a decrease when the strike 
price is low while negative p has the opposite effect. 



20 



Propagators for different p when oc = 

Using 128 time steps and 500.000 configurations 
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Propagators for different p when oc = 0.5 

Using 128 time steps and 500,000 configurations 



3 - 



2 - 



-0.5 




0.5 



Propagators for different p when oc = 1 

Using 128 time steps and 500.000 configurations 
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4.2 The Effect of Mean Reversion on the Option Price 



Several authors including Heston and Hull and White |l5j have considered mean- 
reverting processes as there is some empirical evidence that the volatility follows a mean- 
reverting process. We note that our process also includes mean reversion since A and \x can 
be adjusted so that the volatility performs a mean-reverted process. We find that the effect 
of mean reversion is straightforward in that it only seems to change the implied volatility 
curve so that it moves closer to the mean value. 



Comparison of Mean-Reverted and Non Mean-Reverted Processes 

with a = 0.5, p = -0.5 and mean reversion strength 2 
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Figure 6: Comparison of a mean-reverting process and a non mean-reverting process 

For example, consider figure || In this case, we used an initial volatility of Vq — 0.0625 
(so that a = 0.25). We set A = 0.125 = 2 x 0.0625 and fj, = -2 so that the volatility 
is performing a mean-reverting process with the mean the same as the initial value. We 
see that we indeed obtain the expected behaviour as compared to the non mean-reverting 
process. The mean reverting process just produces an implied volatility curve of a similar 
shape which is closer to the mean value. 

4.3 Calibration with Market Data 

We compare our model with market data to see how well it works. Since volatility information 
is available only as 10 day or 50 day averages and the options we were comparing the market 
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Implied Volatility Curve for European Call Options on the S&P 500 

on Jan 5, 1 998 (maturity : Feb 21 , 1 998) at 3:00 p.m. 
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Figure 7: The implied volatility curve for European call options on the S&P 500 Index 
maturing on Feb. 21, 1998 on Jan. 5, 1998 at 3 p.m. 

data to had 47 days to expiration, we used the initial volatility as a free parameter. The 
market data used were the prices of European call option on the Standard and Poor's 500 
Index at 3 p.m. on Jan 5, 1998 with the maturity date given as Feb 21, 1998. The prices 
taken were either the trade nearest to 3 p.m. if it was within half an hour and the average 
of the bid and ask prices closest to 3 p.m. otherwise. The data are presented in table [l| 
We note that the number of option prices we have is much larger than the number of free 
parameters in the model. The value of the Standard and Poor's 500 Index at the same time 
was given as 965.61 (since the S&P 500 is traded several times every minute on the average, 
obtaining data for it presented no problem). The risk-free interest rate r was 5.131% and 
the annualized dividend yield was 1.617% at the same time. The implied volatility curve for 
the market data is shown in figure |7|. 

Looking at the market data, we immediately see that the implied volatility is almost 
monotonically declining. Hence, according to our model, the correlation is very probably 
negative. We also see that the implied volatilities vary within quite a wide range implying 
that £ must be quite high as the volatility must vary widely for the implied volatility to do 
so. Further, the value of the initial volatility can be seen to be about 25% (according to 
our model, not the actual initial volatility which we could not determine with reasonable 
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Table 1: The prices of European call options on the S&P 500 Index whose maturity was on 
21 Feb, 1998 on Jan 5, 1998 at 3:00 p.m. The prices were taken to be those of the closest 
trade if there was a trade within half an hour and the average of the bid and ask prices 
otherwise. 

accuracy) . 

Since there is no simple functional form for the option price, the calibration was performed 
manually. The reasoning above enabled us to start with fairly accurate values. Thus, while 
we cannot guarantee that the result is the best fit curve in any precise sense, we can see 
from figure |8| that the fit is very good. While the in the money options seem to not fit so 
well, this might be because these options are thinly traded. 

The skeptical reader might comment that the number of free parameters (4) is quite large 

and that a fairly wide range of empirical curves might be fitted. However, we note that our 

model predicts only three possibilities for the implied volatility curve, namely a "smile" (low 

correlation), monotonically increasing (positive correlation) or monotonically decreasing^ 

(negative correlation). The existence of a "frown" or kinks in the implied volatility curve 

would be disastrous for the model. (There does appear to be a small kink in the empirical 

The simulations for non-zero correlation show more interesting behaviour but we are interested only in 
the behaviour for strike prices close to the underlying security price as these are the only kind of options 
traded in the market 
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curve but the scale of the kink is very small and occurs for only one value.) 

Indeed, as shown in figure some calibrated values of the parameters for European call 
otpions on the S&P 500 Index on 21 Feb. 1998 are a = 0.27, p = -0.99, f = 3.7 and a = I. 

Calibration With Market Data 

Obtained Parameters : a = 1 , a = 0.27, p = -0.99 and % = 3.7 

0.29 
0.27 
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"3 
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0.19 
0.17 
0.15 

900 920 940 960 980 1000 1020 1040 
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Figure 8: The above graph shows the implied volatility curve produced by the fitted values 
together with the market data. 

5 Numerical Algorithm for a = 1 

In our previous simulation, we found that the value of a for many typical market data 
is approximately one, that is a ~ 1. It turns out that for this case we have very efficient 
algorithms, namely the bisection method. This algorithm is of intrinsic interest and allows us 
to examine the behavior of stochastic volatility in great detail. In particular, we investigate 
the effect of p on implied volatility and show how smiles turn to frowns and so forth. We 
have intentionally chosen the maturity time to be one year so as to magnify the effects. To 
consider the numerical simulation of our model for a = 1 and A = 0, we first set a = 1 in 
eq© to eq@3|) to get 

< x,y\e- rH \x' >= DY (48) 
J ^(l-p^El^ 




25 



where 



So 
Si 



N X 



e 2 - 



i=l 



(49a) 



if ? 1 

jv ~ 1 * - x> + e B r ~ nf* 

,l-P 2 )eE i= i^ I 7^ 2 



N 

P \ ^ i!i 

7L e2 



£ 2 ' 

Syi + e(/i - —) 



(49b) 
(49c) 



% = Vi-Vi-i, Vn = V = initial volatility. 

In eq(|IS|), we have used the notation J DY to mean product dy I - 

Moreover, the complete information regarding the dynamics and evolution of stock price 
S(t) and its volatility V(t), their cross correlators as well as the fluctuations is given by the 
discrete time path integral equation 



< x, y\e tH \x' > 



dy'p(x,y,r\x',y') 



lim / DXDYe s ' 

AT-»oo 



(50) 



where, for e 
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^(l-p 2 ) 7^ ^2^(1 -p 2 ) 
and where the 'action' S' is defined by 
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-2\2 



5' 



:i-p 2 



V r e yi — -e 2 ( htf c 

e 2 £ v e P 2 s; 



0(e) 



(51) 



To determine numerically the probability P(x,y,r\x' ,y'), one could use the Metropolis 
method || to evaluate numerically the path integral given in eq(|48|) by finding the expec- 

tation value of — , the functional average being performed over configu- 

V / 27re(l-p 2 )Eiie^ 

rations of y n with a probability distribution given by e 5 ° , where So and Si are given by the 
formulae in eq(|4"9|a) and eq(f4*9|b) respectively. However, given the special form of So which 
can be interpreted as the kinetic energy of a free quantum particle, a more efficient method 
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is to generate configurations based on the bisection method |[27|| . In the bisection method, 
unlike the Metropolis case, all configurations generated are accepted and this is therefore a 
better algorithm. 

In the bisection method, the interval r is first divided into half; the sample value of y 

T 1 

at the center is generated by vi^) — ~^{y{fy + 2/( r ) + r ^ 2;Z }? where z is a standard normal 
variate, N(0, 1), namely a normal random variable with mean zero and variance of unity. The 
sub-intervals [0, r/2] and [r/2,r] are further bisected and the above algorithm is repeated 
and so forth. In general, the interval [a, b] is bisected and the sample value at the mid-point 
is given by y( ^ ) = ~{y(a) + y(b) + (b — a) l ^ 2 z], where y(a) and y{b) are the values of 
y(t) at the points a and b respectively. After N bisections, we obtain a time lattice with 2 N 

pSi 

discrete points with spacing e — -^\ and the average of — = is taken over 

^6(1 - P 2)£f =1 e^ 

the configurations generated. 



5.1 Pseudo-Codes for the Algorithm 

Essentially, the main formula that we used to simulate the derivative pricing is eq(^). To 
compute the derivative prices, the program in our simulation can be broken into three major 
steps: 

• Generate 2 N yi variables using the bisection method, 

• Evaluate the path integration using eq(|48D, 

• Compute the derivative price. 

The first step in the computation requires the generation of 2 N variables using the bisec- 
tion method. The pseudo-code for the bisection subroutine is as follows: 

• Initialize an array of size 2^ + 1 (Call it Rand Y(i), for random i/j) 

• Set y\ and yu where M = 2 N + 1 to the initial and end values of the volatility. 

• Do J=l to m 

• Do 1=1 to 2 J ~ X 

• Call Normal Variate iV(0, 1), Norm 
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. Y( M(2'-l) + (2 J -2' + 1 ) ^ , Rand M(/-l) + (2^-l + l) 

r/Vf -i- (2 J ~ 1 — 1 1 / t 

+ Rand Y( ^ >) + Norm(^^)} 

• End of I Loop. 

• End of J Loop. 

Computation of the path integration using eq(pE8|) is straightforward. We first fix x and 
y to be some initial value of stock price and volatility. We then divide the variable x' 
into N divisions^ and store them as arrays. For each value of x' and y', one then calls 
the subroutine Bisect which essentially provides a normalized array of 2 M + 1 points and 
computes using eq([49]b), the average of — . This procedure effectively 

yields the probability in eq(|48|). As a final step, one computes the simulated derivative 
pricing by integrating over y' using a simple integration algorithm based on Bode's rule. 
This integration is exact for any polynomial up to and including degree 5. 



6 Results of Numerical Simulations for a = 1 and A = 

6.1 Parameters 

There are a number of parameters which we can fix freely without affecting the analysis of 
the behaviors in the implied volatility curves. Based on the current market data, we have set 
the interest rate arbitrarily at a plausible rate of 6 % per annum throughout the simulation. 
The other free parameters, which are fixed throughout the simulations, assume the following 
values: 

Interest Rate, r: 0.06 or 6 % 

Time interval, r: 1 year 

Number of bisections: 32 

Number of points for x' integration: 100 

Initial Volatility, V: 1 (i.e. y^ = 0) 



6 Typically, due to the exponential relation of the variables with the stock price, a range between ±7 
should be sufficient. 
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6.2 Simulations 



Under the path integral formalism, one can vary the new parameters like p, £ and p in the 
model and study its effect on the option price. With sufficient simulation, one hopes to 
identify precisely the nature of these changes and interpret their importance to the fluctu- 
ation of the option prices in the market. Like Black-Scholes model, one can also get closed 
form expressions of the option prices under certain constraints and limits. Unlike Black- 
Scholes model, it does not assume constant volatility and thus provides greater flexibility for 
a calibration of the market data. 

We shall see that the variation of the call option prices against strike prices for different 
values of correlation parameter, p can differ significantly from Black-Scholes case. We sim- 
ulate option prices with stochastic volatility using a stock price of 100. In figure |], we have 
plotted the two graphs of call option prices against strike price for a fixed p = 0.001 but 
differing values of £. In this simulation, we have set p = 0. In one graph, we have kept £ as 
0.01 while in the other graph we have fixed the same parameter as 0.001. 

graphicl.gif 

Figure 9: Graph of the call option price with stochastic volatility against strike prices for 
p = and p = 0.001 is plotted. A similar graph computed from Black-Scholes equation is 
also drawn for comparison. 

A graph of the option prices computed from Black-Scholes model with a variance of unity 
is also provided. We observe that despite the small value of p = 0.001, there are still some 
drastic differences between the option prices simulated with stochastic volatility and the 
standard Black-Scholes equation. The difference seems drastic for strike prices below the 
at-the-money position of 100 units but for strike prices above the at-the-money position, 
this difference may not be large. We have checked numerically that the option price with 
stochastic volatility converges, as expected, to that predicted by the Black-Scholes formula 
in the limit £ — > 0. 

We simulate and calculate the implied volatilities for various values of p holding stock 
price constant at 100 and then allowing the strike prices to vary. Our initial simulation 
are based on values of p between 0.1 and 0.9. We believe that one can glean important 
information on the behaviour of the curves with these values. The values of p and £ have 
again been arbitrarily fixed at 0.1. The simulated data for p = 0.1 to p = 0.3 is tabulated 
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and shown in table 0. The implied volatility curves for p = 0.1,0.2,0.3 are plotted against 
different strike prices in figure [It]. The graphs clearly show that the volatility curves assume 
frowns for these positive values of p. 

graphic2.gif 

Figure 10: Graphs of implied volatility curves for positive p values between 0.1 and 0.3. The 
parameters £ and p are set to 0.1 and the stock price is fixed at 100. These curves appear 
as 'frowns'. 
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p = -0.01, £ = 0.01 


Strike Price 


Option 


Volatility 


Option 


Volatility 


Option 


Volatility 


75 


31.2748 


0.635 


38.897 


0.3216 


76.615 


2.16 


80 


30.6564 


0.702 


38.216 


0.443 


74.989 


2.119 


85 


30.0967 


0.749 


37.351 


0.5212 


73.386 


2.08 


90 


29.5342 


0.79 


36.584 


0.5804 


71.726 


2.041 


95 


29.0275 


0.827 


35.907 


0.6298 


69.893 


1.995 


100 


28.5038 


0.861 


35.318 


0.6707 


68.412 


1.967 


105 


27.9804 


0.888 


34.648 


0.7057 


66.718 


1.93 


110 


27.5226 


0.911 


33.992 


0.7377 


65.213 


1.902 


115 


27.0136 


0.931 


33.283 


0.7645 


63.628 


1.871 


120 


26.5956 


0.95 


32.666 


0.7906 


62.134 


1.844 


125 


26.174 


0.971 


32.217 


0.8138 


60.737 


1.821 



Table 2: Simulated data values of call option prices and their implied volatilities for different 
strike prices with varying p values and for a stock price of 100 

The parameter p measures the amount of correlation between the stock price and its 
volatility. Now, correlation coefficient can assume negative values and we should not disre- 
gard this possibility. Also, negative values of p tends to 'whip' up the value of the integrand 
in eq (|49|) or eq(^). We have seen from figure |H] that the volatility curves for positive p 
concave downwards as 'frowns'. At this stage, we may be tempted to surmise that concav- 
ity and convexity of volatility curves are associated with positive and negative values of p 
respectively. We shall soon see that this is not true. 
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We investigate the volatility curves for negative p. For the case in which the stock price 
is fixed at 100 and the parameters p = £ are set to 0.1, the graphs of the implied volatilities 
against strike prices for p = —0.1, —0.2, —0.3 are given in figure [11]. The volatility graph 
for p = —0.1 certainly concave upwards as 'smiles'. However, the graphs for p = —0.2 and 
p = —0.3 flip over and appear as frowns again. Thus, negative values of p do not necessarily 
yield volatility curves which concave upwards as smiles. 

graphic3.gif 

Figure 11: Implied volatility curves for negative p values. The parameters p and £ are fixed 
at 0.1. 

Under path integral formalism, it is always possible to expand the probability given in 
eq(f48[) perturbatively in powers of £ when the values of £ are in general very small [] and 
evaluate exact closed- form expression for the coefficients in the expansion [[UJ. However, 
when £ assumes typically the order of unity, (£ ~ 1 which is what market data indicates) 
it is generally theoretically impossible to perform such an expansion and the perturbative 
method fails. 

However, unlike perturbative analysis, the path integral formalism still permits numerical 
evaluation for these values of £. Some simulated results for large £ = 1 with p — 0.1 are 



shown in figure [L2|. We also investigate the variation of the implied volatility curves against 



strike prices at different initial stock prices. When we fix the values of p and £ to 0.1. with 



p = —0.1, the volatility curves concave upwards as smiles for each stock price. In figure [L3 
we show our results for four different values of stock prices, namely 50, 75, 105 and 200. 
Note that the volatility curves for the last two stock prices, namely 105 and 150, appear to 
coincide with each other. We also note that the large value of p which we have chosen for 



the model may have unnecessarily distorted the graphs in figure [13] and shifted them away 
from the at-the-money position. 

graphic4.gif 

Figure 12: Simulation of Implied Volatility curves for large values of = 1) in which 
perturbation analysis fails. The graphs in this figure are simulated p — 0.1 with varying p 
values. 



7 Strictly speaking, the values of the strike prices and stock prices must also be appropriately tuned before 
a perturbative approach is possible. 
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Figure 13: Graphs of implied volatility against strike price for varying values of stock price, 
S. In this simulation, the parameters p, £ and p have been arbitrarily fixed at 0.1, 0.1, -0.1 

In most mature market, the value of p can be small for some securities. We simulate 
the volatility curves with p set to zero and with the parameters p and £ held small at -0.02 
and 0.01 respectively. In this case, we note two interesting observations. Firstly, negative 
p values do not necessary lead to volatility curves which concave upwards as smiles. As 
seen in figure [15], the concavity of the volatility curves can vary with different stock prices. 
Secondly, if we plot the option prices generated in the simulation with a stock price of 100 



and compare the graph with Black-Scholes model as in figure [14], we observe that the two 
graphs intersect near the strike price of 100. This result is consistent with the observations 
in some models [?] that stochastic volatility appear significant away from the at-the-money 
position. This means that implied volatility is lowest at-the-money and increases as the 
strike prices moves away from that position. 

graphic6.gif 

Figure 14: Graph of option price with p, £ and p set at 0, 0.01, -0.02 and the stock price 
fixed at 100 compared to the Black-Scholes model with the same initial volatility. We have 
also indicated the at-the-money position using a vertical line at the strike price of 100. 

The initial volatility in our simulation has been fixed at 1 so that we can see the full 
effects of volatility to the option prices. Generally we would anticipate a smaller value of 
initial volatility for the market data. We have also performed simulations using a smaller 
initial volatility of 0.2. Figure [16] shows the different implied volatility curves which we have 
obtained using this smaller initial volatility for p = 0.01 and £ = —0.01 with p — at two 
different stock prices, namely S = 50(1) and S = 150(7). On the same graph, we have also 
shown the implied volatility curves for the two stock prices, S = 50(h) and S = 150(h) for 
an initial volatility of unity. We note that the curves can differ significantly for different 
values of initial volatility. Indeed, when we increase the stock price from 50 to 150, we note 
that the implied volatilities generally increase for different strike price at the lower initial 
volatility. However, the reverse effect seems to be taking place when the stock price increase 
over the same range for a high initial volatility. 

Since the results can differ drastically for different initial volatilities, we have also done 
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graphic7.gif 

Figure 15: Volatility curves for different stock price, S. In this simulation, the parameters 
/i, £ and p have been arbitrarily fixed at 0, 0.01, -0.02. 

graphic8.gif 

Figure 16: Volatility curves for two stock prices, S = 50 and S = 150 at different initial 
volatilities. For the lower initial volatility of 0.2, we have denoted the curves by S = 50(1) 
and S = 150(7), whereas for the larger initial volatility of 1, we have labeled the curves 
S = 50(h) and S= 150(h). 



some simulation of the implied volatility curves at an initial volatility of 0.2 for different 
stock prices. Figure [18] and O shows the results of our simulation of the implied volatility 



curves for £ and p fixed at 0.01 and respectively. In figure IB, we have considered a positive 
value of p = 0.01 whereas in figure [T^, we perform the simulation with a negative value of 
p = -0.03. 

graphic9.gif 

Figure 17: Volatility curves for different stock prices between S = 50 and S = 150 with the 
initial volatility set at 0.2. The parameters £, p and p have been set at 0.01, and —0.03 
respectively. 

The concavity of the implied volatility curves can significantly affect the decisions of 
market analysts. In figure O and |18j, one quickly observes the drastic change in the concavity 
as the stock prices vary from 50 to 150. In particular, it is interesting to investigate how 
the implied volatility curves change their concavity between the stock price of 60 and 75 in 



figure [15]. Indeed, figure [L5| shows how a smile can turn into a frown if the initial stock price 
falls from 200 to say 50. 
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Figure 18: Volatility curves for different stock prices between S = 50 and S = 150 with 
the initial volatility set at 0.2. The parameters £, /i and p have been set at 0.01, and 0.01 
respectively. 

graphicll.gif 

Figure 19: Volatility curves for different initial stock prices between S = 63.1 and S = 64.8 
with the parameters fixed at the values for the graph in figure 9. The parameters £, p and 
p have been set at 0.01,0 and —0.02 respectively. 

7 Efficiency of Algorithms 

The normal way of doing Monte-Carlo simulations to find option prices with stochastic 
volatility involves directly simulating the process 

dS = rSdt + WSWdt (52) 
dV = (A + pV)dt + £V a Qdt (53) 

by discretising it using the Euler method to 

AVi = (A + fjV^At + {VfexVAt (54) 
ASi = rSiAt + y/ViSie 2 VAi (55) 

where the standard normally distributed random variables E\ and e 2 with correlation p are 
generated by first generating two uncorrelated standard normally distributed random vari- 
ables 5\ and 62 and using ei = 5i and 62 = p5\ + a/1 — p 2 S 2 - The final values of 5* are stored 
and the option price with strike price K is estimated by considering E[ma.x(S — K, 0)]. This 
is the algorithm used along with the control variate method in Johnson and Shanno [TI| (we 
henceforth call the above algorithm the standard algorithm) . 

Before we compare the efficiency of our algorithm with the standard one, we look at 
the possible sources of error and how the error due to each source scales with run time for 
both. The major source of error for both algorithms is the Monte Carlo error which goes 
as N~ 1 / 2 and which goes as the square root of the run time. Another common source of 
error in both the algorithms is the discretization of the process for y or V. The error in 
y is then of the order of h 1 ^ 2 where h is the time step used. However, the effect of this 
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error on the option price is virtually impossible to estimate and we verify that the number 
of time steps is adequate empirically by comparing two simulations with different numbers 
of time steps. The standard algorithm has a further error due to the discretization of the 
S process. The error is again of the order of h 1 ^ 2 . The effect of this error on the option 
price is of the same order as the payoff of the option is piecewise linear with respect to S. 
Hence, this error also goes inversely as the square root of the run time. Our algorithm has 
other errors due to the finiteness of the integration over x, the interpolation error and the 
quadrature error which goes as h 4 . The error due to the finiteness of the quadrature domain 
can be made completely negligible with minimal effort as the propagator goes as e~ d2 where 
d — x — x'. Hence, for large enough d, this error goes as e~* 2 where t is the run time. This 
is a truly negligible error. The interpolation error is difficult to estimate but we empirically 
show that it is very small for interpolation over 100 points as the propagator is very smooth. 
The quadrature error goes as h 4 (Simpson's rule). Since most of the computer time is spent 
on the Monte Carlo simulation, we can also make this error very small with only a small 
increase in run time. Hence, we see that one of the main advantages of our algorithm is that 
it has a significantly lower error for the same number of configurations (since it has no error 
due to the S simulation since these degrees of freedom have been integrated out). 

When generating a single set of option prices with the same error tolerance, our algorithm 
is about 30 times faster. However, our algorithm has an important advantage in that ti and 
ti are independent of p. Hence, when we generate sets of data with all parameters except 
p fixed, we only need to calculate the new propagator using the terms and integrate over 
the final payoff. This effectively results in an increase in efficiency of a factor of six to seven 
when we calculate the prices for 10 different values of p. Hence, when generating data with 
several values of p, our algorithm is about two hundred times faster. 

At a — 1, we can use an alternative algorithm called the bisection method. Compared to 
usual Monte Carlo methods, the bisection method provides a reasonably fast algorithm for 
the simulation of derivative pricings. We need to integrate the path integral over x' values. 
We first divide the interval for the x' values into 100 divisions. For each value of x', we apply 
the bisection method and divide the interval y into 2 5 = 32 points. This application of the 
bisection method yields one configuration of points. In the algorithm, we sample iV such 
configurations for each value of x'. Using Bode's rule as an integration algorithm, we finally 
integrate over the x' values. In general, we find that it is sufficient to execute only N = 100 
configurations for each value of x' in the program since the gain in numerical precision is not 
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substantial. These results were computed with 100 values of strike prices. 

8 Conclusion 

The path integral formalism provides an additional tool for analyzing derivative pricing. 
Although our model may seem heuristic?], our techniques and computational method is gen- 
erally quite novel in the financial market and it has several added advantages. Firstly, it 
does not assume constant volatility and should be a better tool for analyzing market data. 
Secondly, the additional parameters clearly allow us to study more carefully the behavior of 
the volatilities in the market. And lastly, the formalism provides a reasonably fast algorithm 
for computing the option prices. 

In our simulation, we relied on market data to extract the plausible values of p,, a, p 
and £. Although the value of p in most of the currently well-established markets is believed 
to hover between -0.5 and 0.5, we did not restrict our p to these values. However, it is 
possible for p to assume higher values in more mature markets. Indeed, we have treated 
our simulation largely as an experimental tool to investigate the qualitative behavior for the 
implied volatility. We believe that this new formalism, with proper calibration, can offer a 
formidable tool for getting a more accurate pricing of options. Further, we studied the time 
period from 90 days to one year and set the initial volatility arbitrarily at unity. 

One of the most important qualitative result from our simulations for 1 and a = 1 
is that at-the-money, the implied volatility is equal to the naive fixed volatility of Black- 
Scholes, and the latter is approximately equal to the historical volatility. Deep-in-the-money 
and deep-out-of-the-money show frowns and smiles for implied volatility depends very much 
on the value of p. Finally, the sharp cross-over for the implied volatility for certain special 
values of initial stock price shows that hedging the option for these values could lead to large 
changes in the value of the portfolio and can consequently leads to greater risk. 

Finally, we emphasize that it is very important to calibrate this path integration model 

with the market data for different securities and extract all the relevant information about 

the range of parameters. Such an approach can be extremely interesting and revealing for 

the traders who would like to have some qualitative ideas regarding the actual behavior of 

derivative pricing in the markets. 

8 As pointed out by one referee, the issue of the financial market completeness is rather subtle. Allowing 
trading in stock and riskless bonds implies incompleteness in the market. 



36 



9 Endnotes 

We are grateful to Lawrence Ma (Man-Drapeau, Singapore) and Choon-Peng Toh (Bank 
of America) for discussions and helpful conversations. We also thank G. Bhanot for their 
helpful comments and advices regarding the computational aspect of our work. 
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Appendix 

A A Quantum Mechanical Formulation of the Problem 

We define the Hamiltonian operator as 



2 J dx \ 2 J dy 2 9a; 2 9x9?/ 

t2 e 2j/(a-l) $2 

2 djf' 

(56) 

The Hamiltonian is non- Hermit ian and accounts for the irreversibility of the stochastic pro- 
cesses in finance. From eq(|TT[), we obtain the Merton-Garman-Schrodinger equation 

^ = (r + H(x,y))f, (57) 
f(x,y,T) = max(e x -K,0) (58) 

which can be formally solved as 

POO 

f(x, y, t) = e~ rr / dx'(x, y \ e~ 6r \ x')f(x', T) 



(59) 

f(x, T) = max(e I - K, 0), r = T - t 

While this looks deceptively simple, no analytic solution has been obtained for this equa- 
tion. The special case a — | was solved using a series method by Hull and White |jTH 
and using elementary probability techniques by Heston [13|. the case a = 1 was solved by 
Baaquie [j]. 



B The Lagrangian for the Problem 

The central quantity whose knowledge is sufficient to solve the problem is the conditional 
probability given by the propagator 

P(x, y, T | x', t) = (x,y\ e"^ T | x') (60) 

which can be conveniently handled in the Lagrangian formulation of quantum mechanics. 
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To determine a Lagrangian for the problem, we discretize time so that there are N time 
steps. The time step is then e = j^. The continuous variables x(t) and y(t) are discretized to 
Xi and yi where < i < N. The operator (x, y \ e~ Hr \ x',y') can then be decomposed to 

/oo 
dy'{x,y\{e-^) N \x',y l ) 
■oo 

oo poo poo poo poo 



I 



/OO POO POO POO 

dy N -r ■ dx 1 dy 1 / dy x ( 61 ) 
-oo J —oo J —oo J —oo 



-oo J —oo J —oo J —oo J —oo 

{x N , y N I e~ He I xjv-i, y N -i) ■■■(x 1 ,y 1 | e~ He \ x , y ) 

where x^ = x, y^ = y, x = x' and y = y' . 

We see that if we can find (x,y \ e~ He \ x',y'), we can find the propagator and hence 
the option price. Therefore, let us look at this quantity more closely. Before we consider 
this quantity for the stochastic volatility case, let us consider the Black-Scholes (constant 
volatility) case as it is simpler and retains the essential features. 

In the Black-Scholes case, we only have one variable x (as y is just a constant). We write 

(x | e- kBSt | x) = N BS (e)e lBse (62) 

where N(e) is a normalization constant. We see that 

N BS (e) = — \= (63) 
o"V Aire 

f 1 (8x a 2 \ . ^ 

where 5x = x — x' . 

For the stochastic volatility case, we have 

N(e)e L * = (x,y\e- A <\x',y" 



/dp X f dp y , . _fj . . . \ I I 

27 J 27^' y ' e \P*>Pv)\Px,Pv\x,y 



(65) 



The Hamiltonian in the phase space basis is given by 

A e V , £2 2y(a-l) 



' (e y 6x\ . _ y 5y\ . (66) 

+ \ Y~ r ~~) tPx+ [ 2 ~ ^ - — \ l Py 
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Hence, we have 

w = L £ L * exp r ( 2 ^ + ^' ( - l/2) ^ + (67 , 

<5x e*\ . f5y x _„ C 2 ' 2?/: " 1 \ . ' 

T + r- ~2/ Pie vT e ^ 2 — J Wy 

We obtain in a straightforward but tedious manner 

e v(l/2-a) 

*^ (68) 



and 



+ ?(W)lT + r -2j(T + ^ + ''- 1 ^) (69) 



r 



2(1 -p 2 ) V e 2 
which can be simplified to 

e -y /fo e l/ peW(3/2-a) / 5y £2 2y(a-l) 

L = — — — + r - — - £ — + \e- y + fi 



2(l-p 2 )Ve 2 £ Ve ' 2 

3 2j/(l-a) £2 e 2?;(a-l)N 2 

+ \< " + // 



(70) 



2^ V e 2 

This Lagrangian is difficult to deal with analytically and hence we will consider ways to 
obtain numerical algorithms for the problem. 

It should be emphasized that the above Lagrangian is strictly only correct in the limit 
N — > oo and includes terms of order 0(e) and greater apart from the above expression. 

C Discretized Version of the Action 

The action is defined as 



= J Ldt. (71) 
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The discretized version of the action is given by S = ^YliLi^i + ^( e ) where Li is the 
Lagrangian at time step i. The propagator can be written in terms of the action as 

POO 

(x, y | e-* T \x')= / dy'(x, y \ e~ Ar \ x> \ y 1 ) (72) 



lim / DXDYe s (73) 



where we define 



r>x= e TT ; 



\/27re(l - p 2 ) f = \ J-oo V27re(l-p 2 ) 
~ r 00 dy^ 1 '^ 



DY = r dy [l[ f 



(again x = x', x^ — x, y — y' and yjy = y). We note that the action is quadratic in x. 
This enables us to integrate over the stock price. 
We define 

Q = / DXe Sx = . I / -.e s * 74 

J ^(l-p 2 ) l\J-oo ^(l-p 2 ) 

which is the integral of the action over the stock price. 

We now find Q. The x-dependent term in the Lagrangian is 

L *« = -^37) (t + r - y - —r~ (t + Ae + " - Hr- J J (75) 



Let 



c l = r- — - "—^ i^-f + Ae"* + /i - 5__ ) ( 76) 



Hence, 

N 



S x = -- 7 -^—^y2e- yi (x i -x i . 1 + ec i ) 2 (77) 
2e(l ~ P 2 ) ^ 



We now change the variables to Z; L defined by 



Xj — Zj — e ^y^ j Cj (78) 

3=1 
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Then, 



S z 



1 N 



e K U-^_i) 2 



(79) 



obtaining 







-to/2 



JV-1 

n 



0) 



V27re(l-p 2 ) ^ J-oo V2vre(l - p 

All the integrations can be performed exactly by a process of induction. The exact proce- 
dure can be found in any textbook on path integration such as Kleinert [20| or Roepstorff [[27 
(It is also treated in Baaquie [[[]). We illustrate the method below. 
The integration over z\ is easily performed. We obtain 

dzie- yi ' 2 



-oo v/27re(l - P 2 ) 
e y 2 /2 
, =exp 



exp 



2e(l -P 2 ) 

1 1 



■ w («2-^i) 2 + e- w (zi-z ) 2 



\z% - Zq) 



(81) 



2e(l - p 2 ) e y ^ + e»» 

The above integration can be repeatedly performed over all the variables Zj to obtain 

e Si 



Q 



2^6(1 - p 2 ) Ef =1 ev 



(82) 



where 
5i 



1 


2e(l 


- P 2 ) Ef=i e* 




1 


2e(l 


- P 2 ) Eii 




1 



A' 



i=i 



(83) 



2e(l-p 2 )Ef=i^ 



x 



N 



i=l 



z Vi pe yd3/2-a) / Sy . 



+ Xe~ Vi + fi 



2„2yi(a-l) 



-i 2 



On taking the limit, iV — > oo, we get 
1 



Si 



2(1 -p 2 V 



w p( e ^)( 3 / 2 -") -e^ )^ 2 -")) pA„ pp pr 

x-x + rr — — r - - —6 - + — C 

2 (3/2 -a)£ £ £ 1 2 S 



(84) 
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The term e^ )( 3 / 2 " a ) arises from the fact that £ dte y ^ 2 ~^ § = f v $ dye y ^ 2 ~^ and 



uj= / e y dt= / Vdt (85) 
Jo Jo 

9= e y{1/2 ~ a) dt= / V 1/2 - a dt (86) 
Jo Jo 

v = r e y^ 2 -^dt = r v z ' 2 - a dt (87) 



JO 



C = / e^-^dt = / V a - l ' 2 dt. 
Jo Jo 

Hence, if we can find the joint probability density functions for u, 9, 77, £ and v = V 3 ^ 2 ~ a (0), 
with V following the stochastic process in ([]), we obtain an analytic solution for the problem 
given by 

roo roo poo roc roo ^Si(u>,d,ri,C,,v) 

(x, y | e~ Hr | x') = / du dOj dr] d( dv ===f(u, 6, rj, (, u) 

Jo Jo Jo Jo Jo v27re(l - p A )uj 

(89) 

where / is the joint probability distribution function. 

Hence, we retain the discrete solution which finally gives us 

(x,y | e' Hr | x') = I BY = (90) 

f 2K€(l - p 2 ) Eli <*» 



DYP(x,y,r\x',y') (91) 



where Si is given in (|83| ) and 



So = ~ £ e 2 - (1 - Q) + Ae"» + p - ^ J (92) 

f^t} dv e y ' {1 - a) \ 



We need to start from the path integral given in eq([f3D to study numerically correlators such 
as < e Xn e Ym >. However, if one is interested solely in the price of the option, one needs 
to determine only P(x, y, r\x', y') and in this case, the simplified discrete-time path integral 
obtained in eq(^) should be used as the starting point for numerical studies. 
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